library(org.Mm.eg.db)
library(AnnotationDbi)
library(data.table)
library(tidyverse)
library(rtracklayer)
library(RColorBrewer)
library(gplots)
library(clusterProfiler)
library(fgsea)
library(grid)
dir.create("PDF_figure/KEGG_and_GSEA_analysis_merged", showWarnings = FALSE)
# lead the miRNA list with peaks associated
mirs.peaks <- readRDS("Datafiles/miRNA-merged-peaks-list-09282019-withIDs.rds")
len <- sapply(mirs.peaks, function(x) length(x))
mirs.peaks <- mirs.peaks[order(-len)]
mirna <- names(mirs.peaks)
Run KEGG analysis with the whole universe as background.
KEGG.list <- list()
mirname <- c()
for (i in 1:20) {
sig.genes <- mirs.peaks[[i]]$target_Entrez_ID
kk <- enrichKEGG(gene = sig.genes, organism = 'mmu')
enriched.number <- dim(kk)[1]
print(paste("Number of KEGG pathway enriched for", mirna[i], ":", enriched.number))
if (enriched.number > 0) {
KEGG.list <- c(KEGG.list, list(as.data.frame(kk)))
d <- dotplot(kk, showCategory = 20) + labs(title = paste("Enriched KEGG pathways for", mirna[i], "targets"))
print(d)
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(kk), paste("KEGG Analysis/", str_remove_all(mirna[i], "/"), "_KEGG.csv", sep = ""))
}
}
## Reading KEGG annotation online:
##
## Reading KEGG annotation online:
## [1] "Number of KEGG pathway enriched for let-7-5p/miR-98-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-15-5p/16-5p/195-5p/322-5p/497-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-29-3p : 11"
## [1] "Number of KEGG pathway enriched for miR-200-3p/429-3p : 60"
## [1] "Number of KEGG pathway enriched for miR-17-5p/20-5p/93-5p/106-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-24-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-194-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-27-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-103-3p/107-3p : 1"
## [1] "Number of KEGG pathway enriched for miR-196-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-23-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-26-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-19-3p : 3"
## [1] "Number of KEGG pathway enriched for miR-25-3p/32-5p/92-3p/363-3p/367-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-30-5p/384-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-484 : 3"
## [1] "Number of KEGG pathway enriched for miR-141-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-130-3p/301-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-22-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-96-5p : 0"
pdf("PDF_figure/KEGG_and_GSEA_analysis_merged/KEGG.pdf",
height = 6,
width = 10)
for (i in 1:20) {
sig.genes <- mirs.peaks[[i]]$target_Entrez_ID
kk <- enrichKEGG(gene = sig.genes, organism = 'mmu')
enriched.number <- dim(kk)[1]
print(paste("Number of KEGG pathway enriched for", mirna[i], ":", enriched.number))
if (enriched.number > 0) {
KEGG.list <- c(KEGG.list, list(as.data.frame(kk)))
d <- dotplot(kk, showCategory = 20) + labs(title = paste("Enriched KEGG pathways for", mirna[i], "targets"))
print(d)
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(kk), paste("KEGG Analysis/", str_remove_all(mirna[i], "/"), "_KEGG.csv", sep = ""))
}
}
## [1] "Number of KEGG pathway enriched for let-7-5p/miR-98-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-15-5p/16-5p/195-5p/322-5p/497-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-29-3p : 11"
## [1] "Number of KEGG pathway enriched for miR-200-3p/429-3p : 60"
## [1] "Number of KEGG pathway enriched for miR-17-5p/20-5p/93-5p/106-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-24-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-194-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-27-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-103-3p/107-3p : 1"
## [1] "Number of KEGG pathway enriched for miR-196-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-23-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-26-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-19-3p : 3"
## [1] "Number of KEGG pathway enriched for miR-25-3p/32-5p/92-3p/363-3p/367-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-30-5p/384-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-484 : 3"
## [1] "Number of KEGG pathway enriched for miR-141-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-130-3p/301-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-22-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-96-5p : 0"
dev.off()
## quartz_off_screen
## 2
names(KEGG.list) <- mirname
saveRDS(KEGG.list, "Datafiles/merged.peak.KEGG.result.rds")
I would like to do GSEA analysis on KEGG pathways for each miRNA's associated peaks.
# load the proteomics datafile
protein_DGE <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/Proteomics data/colon tumor-enema model/crcMS_diff.csv")
## Warning: Missing column names filled in: 'X1' [1]
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## X1 = col_double(),
## `Protein Id` = col_character(),
## Description = col_character(),
## p_values = col_double(),
## q_values = col_double(),
## foldChange = col_double(),
## LFC = col_double()
## )
protein_quant <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/Proteomics data/colon tumor-enema model/2017-03-19_Haigis_5v5_Pro.csv")
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## `Protein Id` = col_character(),
## Description = col_character(),
## control1 = col_double(),
## control2 = col_double(),
## control3 = col_double(),
## control4 = col_double(),
## control5 = col_double(),
## KRAS1 = col_double(),
## KRAS2 = col_double(),
## KRAS3 = col_double(),
## KRAS4 = col_double(),
## KRAS5 = col_double()
## )
# annotate the proteomics dataset with Entrez ID
# Entrez IDs are annotated using `org.Mm.eg.db` package since this is the most updated
annotations_entrez <- AnnotationDbi::select(org.Mm.eg.db,
keys = as.character(protein_DGE$`Protein Id`),
columns = c("ENTREZID"),
keytype = "UNIPROT")
## 'select()' returned 1:many mapping between keys and columns
# Determine the indices for the non-duplicated genes
non_duplicates_idx <- which(duplicated(annotations_entrez$ENTREZID) == FALSE)
# Return only the non-duplicated genes using indices
annotations_entrez <- annotations_entrez[non_duplicates_idx, ]
annotations_entrez <- annotations_entrez[!is.na(annotations_entrez$ENTREZID),]
# Check number of NAs returned
is.na(annotations_entrez$ENTREZID) %>%
which() %>%
length()
## [1] 0
# annotate the dataset with Entrez ID
protein_DGE$Entrez_ID <- NA
for (i in 1:dim(protein_DGE)[1]) {
index <- grep(protein_DGE$`Protein Id`[i], annotations_entrez$UNIPROT)
if (length(index)>0) {
protein_DGE$Entrez_ID[i] <- paste(annotations_entrez$ENTREZID[index], collapse = " ")
}
}
# calculate signal-to-noise ratio for GSEA later
s2n <- function(num_list, cond_1 = c(4:6), cond_2 = c(1:3)) {
mean1 <- mean(as.numeric(num_list[cond_1]))
if (mean1 == 0) {
mean1 = 1
}
mean2 <- mean(as.numeric(num_list[cond_2]))
if (mean2 == 0) {
mean2 = 1
}
sd1 <- sd(as.numeric(num_list[cond_1]))
sd2 <- sd(as.numeric(num_list[cond_2]))
sd1 <- min(sd1, 0.2*abs(mean1))
sd2 <- min(sd2, 0.2*abs(mean2))
s2nvalue <- (mean1-mean2)/(sd1+sd2)
return(s2nvalue)
}
protein_DGE$s2n <- apply(protein_quant,1,s2n,cond_1 = c(8:12),cond_2 = c(3:7))
# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_H_v5p2.rdata")
pathways <- Mm.H
# GSEA on Hallmark gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
gseadata <- mirs.peaks[[i]]$target_Entrez_ID
gseadata <- gseadata[!is.na(gseadata)]
ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
if (dim(fgseaRes)[1] > 0) {
GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_Hallmark_mergedpeaks.csv", sep = ""))
for (n in 1:length(fgseaRes$pathway)) {
p <- plotEnrichment(pathways[[fgseaRes$pathway[n]]], ranks) + labs(title = paste(mirna[i],"-",fgseaRes$pathway[n], sep = ""))
grid.newpage()
print(p)
}
grid.newpage()
plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5)
}
}
pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_Hallmark_merged_peaks.pdf',
width = 10,
height = 6)
for (i in 1:20) {
gseadata <- mirs.peaks[[i]]$target_Entrez_ID
gseadata <- gseadata[!is.na(gseadata)]
ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
if (dim(fgseaRes)[1] > 0) {
GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_Hallmark_mergedpeaks.csv", sep = ""))
for (n in 1:length(fgseaRes$pathway)) {
p <- plotEnrichment(pathways[[fgseaRes$pathway[n]]], ranks) + labs(title = paste(mirna[i],"-",fgseaRes$pathway[n], sep = ""))
grid.newpage()
print(p)
}
grid.newpage()
plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5)
}
}
dev.off()
## quartz_off_screen
## 2
names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.Hallmark.result.rds")
# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_c2_v5p2.rdata")
pathways <- Mm.c2
# GSEA on C2 gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
gseadata <- mirs.peaks[[i]]$target_Entrez_ID
gseadata <- gseadata[!is.na(gseadata)]
ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
if (dim(fgseaRes)[1] > 0) {
GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C2_mergedpeaks.csv", sep = ""))
topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
print(mirna[i])
grid.newpage()
plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5)
}
}
## [1] "let-7-5p/miR-98-5p"
## [1] "miR-15-5p/16-5p/195-5p/322-5p/497-5p"
## [1] "miR-29-3p"
## [1] "miR-200-3p/429-3p"
## [1] "miR-196-5p"
## [1] "miR-26-5p"
## [1] "miR-19-3p"
## [1] "miR-30-5p/384-5p"
## [1] "miR-484"
pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_C2_merged_peaks.pdf',
width = 10,
height = 6)
for (i in 1:20) {
gseadata <- mirs.peaks[[i]]$target_Entrez_ID
gseadata <- gseadata[!is.na(gseadata)]
ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
if (dim(fgseaRes)[1] > 0) {
GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C2_mergedpeaks.csv", sep = ""))
topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
print(mirna[i])
grid.newpage()
plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5)
}
}
## [1] "let-7-5p/miR-98-5p"
## [1] "miR-15-5p/16-5p/195-5p/322-5p/497-5p"
## [1] "miR-29-3p"
## [1] "miR-200-3p/429-3p"
## [1] "miR-196-5p"
## [1] "miR-26-5p"
## [1] "miR-19-3p"
## [1] "miR-30-5p/384-5p"
## [1] "miR-484"
dev.off()
## quartz_off_screen
## 2
names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.C2.result.rds")
# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_c3_v5p2.rdata")
pathways <- Mm.c3
# GSEA on C3 gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
gseadata <- mirs.peaks[[i]]$target_Entrez_ID
gseadata <- gseadata[!is.na(gseadata)]
ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
if (dim(fgseaRes)[1] > 0) {
GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
mirname <- c(mirname,mirna[i])
write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C3_mergedpeaks.csv", sep = ""))
topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
print(mirna[i])
grid.newpage()
plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5)
}
}
## [1] "let-7-5p/miR-98-5p"
## [1] "miR-29-3p"
## [1] "miR-200-3p/429-3p"
## [1] "miR-17-5p/20-5p/93-5p/106-5p"
## [1] "miR-24-3p"
## [1] "miR-27-3p"
## [1] "miR-103-3p/107-3p"
## [1] "miR-196-5p"
## [1] "miR-19-3p"
## [1] "miR-141-3p"
## [1] "miR-22-3p"
names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.C3.result.rds")
# compile a matrix containing all genesets that had enrichment for all 20 miRNAs and their enrichment score
hall.GSEA <- readRDS("Datafiles/merged.peak.GSEA.Hallmark.result.rds")
c2.GSEA <- readRDS("Datafiles/merged.peak.GSEA.C2.result.rds")
#c3.GSEA <- readRDS("Datafiles/merged.peak.GSEA.C3.result.rds")
gsea.matrix <- data.frame(c(1:20))
gsea.matrix <- t(gsea.matrix)
colnames(gsea.matrix) <- mirna[1:20]
gsea.matrix[1,] <- NA
for (i in 1:length(hall.GSEA)) {
for (n in 1:length(hall.GSEA[[i]]$pathway)) {
if (length(grep(hall.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
col.index <- grep(names(hall.GSEA[i]), colnames(gsea.matrix))
row.index <- grep(hall.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
gsea.matrix[row.index, col.index] <- hall.GSEA[[i]]$ES[n]
}
else {
gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- hall.GSEA[[i]]$pathway[n]
col.index <- grep(names(hall.GSEA[i]), colnames(gsea.matrix))
gsea.matrix[dim(gsea.matrix)[1], col.index] <- hall.GSEA[[i]]$ES[n]
}
}
}
gsea.matrix <- gsea.matrix[-1,]
for (i in 1:length(c2.GSEA)) {
for (n in 1:length(c2.GSEA[[i]]$pathway)) {
if (length(grep(c2.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
col.index <- grep(names(c2.GSEA[i]), colnames(gsea.matrix))
row.index <- grep(c2.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
gsea.matrix[row.index, col.index] <- c2.GSEA[[i]]$ES[n]
}
else {
gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- c2.GSEA[[i]]$pathway[n]
col.index <- grep(names(c2.GSEA[i]), colnames(gsea.matrix))
gsea.matrix[dim(gsea.matrix)[1], col.index] <- c2.GSEA[[i]]$ES[n]
}
}
}
#for (i in 1:length(c3.GSEA)) {
# for (n in 1:length(c3.GSEA[[i]]$pathway)) {
# if (length(grep(c3.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
# col.index <- grep(names(c3.GSEA[i]), colnames(gsea.matrix))
# row.index <- grep(c3.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
# gsea.matrix[row.index, col.index] <- c3.GSEA[[i]]$ES[n]
# }
# else {
# gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
# rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- c3.GSEA[[i]]$pathway[n]
# col.index <- grep(names(c3.GSEA[i]), colnames(gsea.matrix))
# gsea.matrix[dim(gsea.matrix)[1], col.index] <- c3.GSEA[[i]]$ES[n]
# }
# }
#}
gsea.matrix <- as.matrix(gsea.matrix)
my_palette <- colorRampPalette(c("blue", "white", "red"))(128)
na_color <- colorRampPalette("grey")
png('GSEA Analysis/GSEA_HEAP-CLIP_merged_peaks.png',
width = 1200,
height = 800,
res = 100,
pointsize = 8)
par(cex.main=1.8)
heatmap.2(gsea.matrix,
main = "GSEA for 20 most active\nmiRNAs from merged HEAP-CLIP",
density.info = "none",
key = TRUE,
lwid = c(1,7),
lhei = c(1,7),
col=my_palette,
labRow = rownames(gsea.matrix),
margins = c(20,25),
symbreaks = TRUE,
trace = "none",
Rowv=FALSE,
Colv=FALSE,
dendrogram = "none",
na.color = "grey",
cexCol = 1,
na.rm = FALSE)
dev.off()
## quartz_off_screen
## 2
pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_HEAP-CLIP_merged_peaks.pdf',
width = 24,
height = 16)
par(cex.main=1.8)
heatmap.2(gsea.matrix,
main = "GSEA for 20 most active\nmiRNAs from merged HEAP-CLIP",
density.info = "none",
key = TRUE,
lwid = c(1,7),
lhei = c(1,7),
col=my_palette,
labRow = rownames(gsea.matrix),
margins = c(20,25),
symbreaks = TRUE,
trace = "none",
Rowv=FALSE,
Colv=FALSE,
dendrogram = "none",
na.color = "grey",
cexCol = 1,
na.rm = FALSE)
dev.off()
## quartz_off_screen
## 2
GSEA